Short-time Rheology and Diffusion in Suspensions of Yukawa-type Colloidal Particles 
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A comprehensive study is presented on the short-time dynamics in suspensions of charged colloidal 
spheres. The explored parameter space covers the major part of the fluid-state regime, with colloid 
concentrations extending up to the freezing transition. The particles are assumed to interact directly 
by a hard-core plus screened Coulomb potential, and indirectly by solvent-mediated hydrodynamic 
interactions (His). By comparison with accurate accelerated Stokesian Dynamics (ASD) simulations 
of the hydrodynamic function H{q), and the high-frequency viscosity rjoo , we investigate the accuracy 
of two fast and easy-to-implement analytical schemes. The first scheme, referred to as the pairwise 
additive (PA) scheme, uses exact two-body hydrodynamic mobility tensors. It is in good agreement 
with the ASD simulations of H{q) and rjoo, for smaller volume fractions up to about 10% and 20%, 
respectively. The second scheme is a hybrid method combining the virtues of the ^7 scheme by 
Beenakker and Mazur with those of the PA scheme. It leads to predictions in good agreement with 
the simulation data, for all considered concentrations, combining thus precision with computational 
efficiency. The hybrid method is used to test the accuracy of a generalized Stokes-Einstein (GSE) 
relation proposed by Kholodenko and Douglas, showing its severe violation in low salinity systems. 
For hard spheres, however, this GSE relation applies decently well. 



I. INTRODUCTION 

Charge-stabilized systems of globular Brownian par- 
ticles are ubiquitously found over a large range of par- 
ticle sizes, from large, micron-sized colloids [1-3] down 
to nanometer-sized proteins [4-6]. For many such sys- 
tems, where van der Waals attractions are to a good ap- 
proximation negligible, the pair interactions can be de- 
scribed to good accuracy by a hard-sphere plus repul- 
sive Yukawa (HSY) type pair potential, of range deter- 
mined by the ionic strength of dissolved co- and coun- 
terions. The HSY model spans the range from neu- 
tral hard spheres, corresponding to zero screening length 
or vanishing Yukawa potential strength, to long-range 
electric repulsion occurring in low-salinity systems with 
large screening length. In most applications of the 
HSY model to charged colloids, the Yukawa-tail of the 
pair interaction is described by the electrostatic part of 
the Derjaguin-Landau-Vcrwey-Overbeck (DLVO) poten- 
tial [7]. 

In experimental data analysis and many theoretical ap- 
plications, easy-to-implement analytic methods are on 
demand that allow to calculate, with good accuracy, 
short-time dynamic properties, such as diffusion coeffi- 
cients and high-frequency viscosities with a minimal com- 
putational effort. Short-time diffusion properties are rou- 
tinely measured in dynamic light scattering [8, 9], X-ray 
photon correlation spectroscopy [2, 10], and neutron spin 
echo [11, 12] experiments. The high-frequency viscosity 
Tjoa is probed experimentally using torsional viscometers, 
or on employing approximate generalized Stokes-Einstein 
(GSE) relations, which relate rjao to a diffusion property 
[13, 14]. 
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Fast and accurate theoretical methods for calculat- 
ing dynamic properties are valuable in particular for an 
extensive data analysis, where different system parame- 
ters such as concentration, salt content, particle size and 
charge, pH-value, and solvent properties, are considered 
in a broad range of values. 

We point out here that short-time dynamic proper- 
ties are of relevance not only in their own right. They 
are also required as input to theories describing colloidal 
long-time dynamics such as mode-coupling and dynamic 
density functional theory approaches. 

A particular challenge in the development of analytic 
methods is the inclusion of the long-ranged, and for 
more concentrated systems non-pairwise additive, hydro- 
dynamic interactions (His), which essentially influence 
the dynamic properties in their short-time behavior. Due 
to their complexity, the inclusion of His constitutes a se- 
vere bottleneck in Brownian dynamics simulations. The 
account of many-body His in analytic methods is only 
possible by introducing approximations. For this rea- 
son, it is of prime importance to assess the overall accu- 
racy of analytic methods, by the comparison to precise 
benchmark results obtained from computationally elab- 
orate dynamic computer simulations. 

In this paper, we discuss the pros and cons of two easy- 
to-implement analytic methods of calculating short-time 
dynamic properties, such as the wavenumber-dependent 
hydrodynamic function H(q), and the low shear-rate, 
high-frequency limiting viscosity 7700- The first method, 
referred to as the pairwise additive (PA) scheme, uses 
exact two-body hydrodynamic mobility tensors includ- 
ing the lubrication part, but it fully disregards three- 
body and higher-order hydrodynamic contributions. The 
second scheme is a hybrid method combing the virtues 
of Beenakker and Mazur's so-called 57-scheme approach 
for H{q) [1-5, IG] and rjoo [IT], with those of the PA 
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scheme, and precise known results for the special case 
of neutral hard spheres. The 5^ scheme accounts for 
many-body HI contributions in an approximate way. We 
present the two methods in a self-contained way, allow- 
ing for their easy implementation. Both methods re- 
quire the static structure factor, S{q), or equivalently 
the radial distribution function (rdf) g{r)^ as the only 
input. We calculate this static input using our recently 
developed analytic modified penetrating-background cor- 
rected rescaled mean spherical approximation (MPB- 
RMSA) scheme [18, 19], which allows for a fast and ac- 
curate evaluation of the HSY pair-structure functions. 

The accuracy of both methods for calculating H{q) 
and rjoo is assessed through comparison with a large num- 
ber of simulation results, representative of the full fluid- 
state regime, which we have obtained using accelerated 
Stokesian Dynamics simulations. The usefulness of the 
^7 scheme based hybrid method is illustrated by testing 
the validity of three generalized Stokes-Einstein relations. 

The paper is organized as follows. Sec. II explains 
the essentials of the HSY model, and the MPB-RMSA 
method of calculating S{q). The theoretical background 
on the short-time dynamics of interacting colloidal parti- 
cles is included in Sec. HI, and Sec. IV explains the em- 
ployed methods of calculating short-time dynamic prop- 
erties. Our results for H{q), rjao and additional related 
short-time properties are summarized in Sec. V. Sec. VI 
includes the test of GSEs, notably that proposed by 
Kholodenko and Douglas. Our conclusions are given in 
Sec. VII. 



II. PAIR-POTENTIAL AND STATIC 
STRUCTURE 

The present study is concerned with charged spheri- 
cal colloidal particles that interact directly via the hard- 
sphere plus Yukawa (HSY) repulsive pair potential 
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The coupling amplitude 7 and the screening parameter 
k are given by 



Lb f Ze^'^ 



1 — 
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This constitutes the repulsive part of the DLVO poten- 
tial [7]. Here, /? = l/ksT, with Boltzmann constant 
fcs, absolute temperature T, colloidal hard-core diame- 
ter cr, solvent-characteristic Bjerrum length Lb = /?e^/e 
in Gaussian units, proton elementary charge e, and sol- 
vent dielectric constant e. In the DLVO expressions in 
Eqs. (2) and (3), 7 and k are determined by the number 



concentration of monovalent (salt) colons, Ug, the effec- 
tive colloid charge, Ze, and the colloid volume fraction, 
(p. The square of A; is a sum of two contributions, namely 
fc2 = 2401 Z I Lb /[ct(1 - (j))] and k^ = Smisa^Ls/il ~ 4>). 
The first one, fc^, describes the screening influence of 
countcrions released from the colloid surfaces, and the 
second one, fc^, accounts for the screening influence of 
monovalent electrolyte ions arising from added salt. As 
shown in Refs. [20, 21], the factor 1/(1 - 0) in Eq. (3) 
corrects for the free volume accessible to the microions 
in presence of impermeable colloidal spheres. 

In recent work [18, 19], we have introduced a semi- 
analytic Ornstein-Zernike integral equation scheme [22] 
for calculating equilibrium pair-distribution functions, 
denoted as the modified penetrating-background cor- 
rected mean spherical approximation (MPB-RMSA). 
This computationally highly effective method allows for 
calculating the static structure factor, S{q), as a func- 
tion of the scattering wavenumber g, for particles inter- 
acting by the HSY pair potential given in Eq. (1). The 
MPB-RMSA is a simple improvement of the PB-RMSA 
method by Snook and Hayter [2.H], which in turn is based 
on the frequently used RMSA method by Hansen and 
Hayter [24]. The analytic simplicity, and the low com- 
putational cost of the standard RMSA method is pre- 
served in the MPB-RMSA. In addition, the MPB-RMSA 
constitutes a significant improvement over the RMSA by 
correcting for the RMSA-typical underestimation of the 
principal peak height of S{q). The excellent accuracy of 
the MPB-RMSA has been assessed in Ref. [IN] by means 
of extensive parameter studies in comparison to Monte- 
Carlo (MC) simulations, and results obtained from the 
highly accurate, but non-analytic Rogers- Young (RY) in- 
tegral equation scheme [25]. The additional virtue of the 
MPB-RMSA, which we take advantage of in the present 
study, is its fast performance. For a given parameter 
set {7, fc, 0}, it allows to compute S{q) in an extended 
q range in about 0.1 seconds of cpu time on a desktop 
PC, which is orders of magnitude faster than using the 
RY-scheme, or even more time consuming simulations. 

In the limiting case of neutral hard spheres (HS), at- 
tained for 7 = (Z = 0) or fc — CO (very large ris), 
the MPB-RMSA reduces to the analytic Percus-Yevick 
solution [20, 27], which is known to give accurate pair- 
correlation functions for < 0.35. For larger values 
0.35 ^ < 0.49, the Percus-Ycvick solution tends to 
overestimate somewhat the structure factor peak value, 
S{qm), of hard spheres, and to underestimate the rdf con- 
tact value g{x = 1+). 

Most of the results on short-time properties discussed 
in this work are for the parameters Lb = 5.617 nm, 
a ~ 200 nm, and Z ~ 100, representative of strongly 
charged colloidal spheres suspended in an organic solvent. 
Numerous MPB-RMSA, RY, and MC simulation results 
for S(q) and g{r) using these parameters are included in 
Figs. 2, 3 and 4 of Ref. [IS]. These results cover basi- 
cally the whole fluid regime, with (j) ranging from 10~^ 
to 0.15, and Us from to 10~^ M. For conciseness, and 
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since these results for S{q) and g{r) have been pubUshed 
aheady, we refrain from repfotting them in the present 
work. 



III. SHORT-TIME DYNAMIC PROPERTIES 

For characterizing the coUoidal short-time regime, 
one considers the momentum relaxation time, tb ~ 
TO/(37r?7ocr), the time scale th = cr^ ps / i'^'Ho) of hydro- 
dynamic vorticity diffusion, and the interaction time 
T/ = o'^/(4do) [28, 29], where m is the mass of a colloidal 
sphere, dg = ksT / {^irrjoa) is the translational free dif- 
fusion coefficient for stick hydrodynamic surface bound- 
aries, and PS and rjo are the mass density and shear vis- 
cosity of the suspending Newtonian solvent, respectively. 
For a coarse-grained time-resolution where t ^ tb th, 
the motion of a colloidal particle is erratic and over- 
damped. In the present work we focus on the colloidal 
short-time regime tb <IC t "C t/, during which a particle 
has moved a tiny fraction of its size only. This allows for 
calculating short-time properties using pure equilibrium 
averages. 

Diffusion properties can be measured by a variety of 
scattering techniques, which commonly determine the 
dynamic structure factor [22], 

Siq, t)=\im/^^J2 • [Ri (0) - (t)] } ^ , (4) 

as a function of scattering wavenumber q and correlation 
time t. The brackets, < ... >, denote an equilibrium en- 
semble average. N is the number of colloid particles in 
the scattering volume, q is the scattering wave vector, 
and Rn(i) is the position vector pointing to the center 
of the n-th colloidal particle at time t. Moreover, limoo 
denotes the thermodynamic limit — > oo and system 
volume V — > oo, with n = N/V fixed, which character- 
izes a macroscopic system. On the colloidal short-time 
scale, S{q, t) decays exponentially according to [s] 



S{q,t) 



exp [-q^D{q)t] , 



(5) 



where D{q) is the wavenumber-dependent short-time dif- 
fusion function. A statistical-mechanical expression for 
D{q) follows from the generalized Smoluchowski equation 
in form of the ratio [28-30] 



D{q) = do 



Hjq) 
S{q)- 



(6) 



of the hydrodynamic function 



N 



H{q) = hm ( ^ ^ q . /.«(R^) • qexp {*q • [R, - R,]} 



and the static structure factor S{q) = S{q^t = 0). Here, 
q is the unit vector in the direction of q, and /.i**(R^) is a 
translational mobility tensor linearly relating the hydro- 
dynamic force on a sphere j to the translational velocity 
of a sphere I. This mobility depends in general on the 
instantaneous positions, R^, of all particles through 
the specified hydrodynamic boundary conditions. In this 
work, stick hydrodynamic boundary conditions are as- 
sumed throughout. The positive- valued hydrodynamic 
function H{g) is a measure of the influence of His on 
short-time diffusion. In the (hypothetical) case of hydro- 
dynamically non-interacting particles, H[q) = 1, inde- 
pendent of q. The hydrodynamic function can be inter- 
preted as the reduced short-time generalized mean sedi- 
mentation velocity measured in a homogeneous suspen- 
sion subject to a weak force field coUinear with q and 
oscillating spatially as cos(q • r). Hence, 



lim H{q) 



Us. 



Uo 



(8) 



is equal to the concentration-dependent (short-time) sed- 
imentation velocity. Used, of a slowly settling suspension 
of spheres in units of the sedimentation velocity, Uq, at 
infinite dilution. 

The function H{q) can be expressed as the sum. 



Hiq) 



ds_ 

do 



H'{q), 



(9) 



of a g-dependent distinct part, H'^{q), which van- 
ishes for g — >■ oo, and a self- part equal to the re- 
duced short-time translational self-diffusion coefficient 
ds/do- The diffusion coefficient dg is equal to the 
short-time slope of the mean-squared displacement, 
W{t) = 1/6 < [R{t) - R(0)]^ >, of a colloidal particle 
[2')]. 

Two additional diffusion coefficients related to D{q) 
are the short-time collective diffusion coefficient dc = 
doK/S{q —7- 0), and the short-time cage diffusion coef- 
ficient dcgc = D{qm)- These two coefficients characterize 
the decay rates, respectively, of thermally induced con- 
centration fiuctuations of macroscopic wavelengths, and 
of a wavelength related to the size, l-Kjqmi of the dynamic 
next-neighbor cage formed around a particle. 

So far only diffusion properties have been discussed. 
A rheological short-time property is the high-frequency 
limiting viscosity, ryoo, which linearly relates the average 
deviatoric suspension shear stress to the applied rate of 
strain in a low-amplitude, oscillatory shear experiment 
with frequency uj 3> 1/t/. The statistical- mechanical 
expression for 7700 is [31] 

'?oo = ^o + lhn-l- ^ ljZl4t^,,S^''^^ (10) 



Q,/3 = l \l,j = \ 



invoking the Cartesian components, a'j'f 



, of the 3 X 



(7) 3x3x3 dipole-dipole mobility tensor /xf", that relates 
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the symmetric hydrodynamic force dipole moment tensor 
of sphere I to the rate of strain tensor related to sphere 
j. For stick hydrodynamic boundary conditions, rjao = 
r]o[l + 2.50 + 0((/)^)], where 770 is the solvent viscosity. 
The 0(0^) contribution is due to particle interactions. 
Note here that < ... > describes an equilibrium average 
with respect to the unsheared system. Direct interactions 
affect r]oo only through their influence on the equilibrium 
particle distribution. 

The great difficulty in evaluating Eqs. (7) and (10) 
to obtain H{q) and 7700, respectively, lies in the calcula- 
tion of the hydrodynamic tensors /x**(R^) and fif^{Il^), 
and in the associated many-particle average. Except for 
numerically expensive simulations [14, 32-34], it is prac- 
tically impossible to gain numerically precise results for 
concentrated systems. In searching for methods which 
allow to calculate H(q) and rj^ to decent accuracy with 
moderate numerical effort, one has to resort to approx- 
imate methods. Two of these methods, namely the so- 
called pairwise additive (PA) approximation, and the ^7- 
scheme by Beenakker and Mazur supplemented by a so- 
called self-part correction, are discussed in the next sec- 
tions. The methods are presented in a self-contained way 
to facilitate their implementation by an interested reader. 
Both methods have in common that they require S{q), or 
equivalently g{r), as the only input. The pros and cons 
of both methods are assessed in comparison to elaborate 
computer simulations. 



IV. COMPUTATIONAL METHODS 



A. Pairwise additive approximation 



In the PA approximation, the TV-particle translational 
mobility tensors, /xJ*(R^), are approximated by the sum 
of two-body mobilities according to 



ksT 
do 



PA 



N 



l+^an(Rz-R„ 
(l-%)ai2(Ri-Rj)- (11) 



The 3x3 mobility tensor 1 -I- an relates, for an isolated 
pair of particles in a quiescent fluid, the force on particle 
1 to its own velocity. Correspondingly, a.12 relates the 
force on particle 2 to the velocity of particle 1. The ax- 
ial symmetry of the two-sphere problem allows to split 
the reduced mobilities into longitudinal and transverse 
components. 



a,,(r) ==x^,(r)ff-K2/^,(r) [1 



(12) 



where we use the notation from [35]. The mobility com- 
ponents x^j{r) and Uijir) can be expanded analytically 
in powers oi a/r = 1/x using recursion formulas [3()]. 

In a homogeneous fluid system, the ensemble average 
of a function / depending on two particle coordinates can 



be expressed in the thermodynamic limit as 

1 



Rj)) = fim 



V 



dr g(r)/(r). 



(13) 



The combination of Eqs. (7), (11), and (13) leads to the 
following PA results for ds, 



ds_ 
do 



1 



PA 



dxx^gix) [xl,{x)+2y1^{x)-3], 

(14) 



and for the distinct part of the hydrodynamic function, 

Jiiy) , 



H'{y)\ 



PA 



450- 



180 / dx xh{x) 



y 

jo{xy) 



(15) 



ji{xy) , j2{xy) 



xy 



240 
240 



dx x^g{x)yl2{x)io{xy) + 



dx x^g{x) [il2{x) - y'l2{x)] X 



h{xv) 



xy 



Here, y = qa is the diameter-scaled wavenumber, j„ is 
the spherical Bessel function of first kind and order n, 
and h is the total correlation function defined by h{x) ~ 
gix)-l. 

We have introduced here the short-range mobility 
parts 



yi2ix)= yi2ix) 



' 3/4a;"i 
3/8x^1 - 



- l/8a;~^ 
l/16x-^ 



(16) 
(17) 



which include all terms in the series expansion except for 
the far-field terms up to the dipolar (Rotne-Prager) level, 
which are subtracted off. 

Analogous to the translational diffusivity tensors, the 
dipole-dipolc mobility tensors are approximated in the 
PA scheme by their self and distinct two-body contribu- 

tions /x'^'^ii M'''*i2 (i"): respectively. In hydrodynami- 
cally semi-dilute suspensions, the high-frequency viscos- 
ity of colloidal spheres at low shear-rate is then obtained 
from [14, 37, 38] 



14. 

-^(^^^ 25^ 



(1 



dx x^g{x)J{x), (18) 



,,<i(i(2) / N 



(19) 



where the rapidly decaying two-body shear mobility func- 
tion, J(a;), accounts for the two-body His. For stick hy- 
drodynamic boundary conditions, J{x) = 15/128 x~^ + 
0(x"8). 

For pair-distances x > 3, we use explicit analytic ex- 
pansions up to 0(a;~20)^ given in Ref. [30] for the two- 
body mobility functions xf^ and y^^ and the leading- 
order far-field expression, J{x) = 15/128 a;~^, for the 
shear mobility function. Since the expansions in 1/x 
converge slowly at small separations, we employ accurate 
numerical tables for a; < 3, which in particular account 
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for lubrication at near-contact distances. The tables are 
based on recursion expressions and a lubrication analysis 
given in [35]. 

Using the zeroth-ordcr concentration-expansion for the 
rdf of hard spheres given by g^^{x) = Q{x — 1) + 0(0), 
with Q denoting the unit step function, we have checked 
that our PA code precisely reproduces the truncated 
virial expressions 



1-1.83150-1-0(02), (20) 
1 - 6.5460-^0(02), (21) 
1 + I0 + 5.002302 -f 0(03), (22) 



which were obtained in [39-41] using a mobility series 
expansion including terms up to 0(a;~^°™), and a lubri- 
cation correction. 

All PA results for H[q) and 7700 discussed in this paper 
are based on Eqs. (14), (15) and (18), with g{x) com- 
puted in MPB-RMSA. Dynamical properties predicted 
by the PA scheme are exact to linear order in 0. Thus, the 
PA scheme is especially well-suited for hydrodynamically, 
but not necessarily structurally, dilute systems. Charge- 
stabilized suspensions at low salinity and concentration, 
where near-contact configurations are very unlikely, are 
prime examples of hydrodynamically dilute, but struc- 
turally ordered systems, showing pronounced oscillations 
in S(ci) and g{r). 

Moreover, the PA scheme can be used to check the 
accuracy of other approximate schemes, such as the 
(self-part corrected) 5^ scheme, in the low concentra- 
tion regime. At larger volume fractions, however, and 
for diffusion properties like dc and if, where non-pairwise 
additive His are particularly influential, the PA approx- 
imation is bound to fail. Note that, while in the present 
work particles with stick hydrodynamic boundary condi- 
tions are considered, the PA scheme can be easily gen- 
eralized to porous particles, and particles with slip-stick 
boundary conditions, simply by using the corresponding 
two-body mobility functions given, e.g., in [36]. 



B. (57-method by Beenakker and Mazur 

Different from the PA scheme, which cannot be applied 
to concentrated systems, the renormalized concentration 
fluctuation (termed ^7) expansion method of Beenakker 
and Mazur [15, 17] is applicable to fluid-ordered systems 
also at large values of 0, where three-body and higher- 
order HI contributions arc important. The ^7 method is 
an effective medium approach based on a partial resum- 
mation of many-body HI contributions. While applicable 
to all 0, its results for H{q) and 7700 reveal moderate inac- 
curacies at all concentrations, including the very dilute 
regime where the PA approach becomes exact. These 
inaccuracies can be partially traced back to the approx- 
imate expressions of /x**(R^) and fifHlt'^) used in the 



derivation of the 57-scheme where, in particular, lubrica- 
tion corrections are disregarded. Higher order terms in 
the ^7 expansion require as input static correlation func- 
tions of increasing order (pair, triplet, and so on) with 
swiftly increasing difficulties in their evaluation. 

In the present study, we use the easy-to-implement 
standard version of the Sj method for which (like for the 
PA scheme) only S{q) is required as input, with the latter 
calculated here using the MPB-RMSA scheme. This is 
the zeroth-order ^7 approximation regarding H{q), and 
the second-order approximation regarding 7700 • 

The zeroth-order ^7 scheme for H{q) has been applied 
in the past both to neutral and charged colloidal parti- 
cles, but the second-order scheme for 7700 was used so 
far for neutral hard spheres only. To our knowledge, the 
present work provides the first test of the scheme for 
charged, Yukawa-type particles. 

The zeroth-order J7-scheme expression for H{q) con- 
sists of a microstructurc- independent part. 



dM 



do 





sm{y) 


57 ^ Jo 


. y . 



[l + 057„(y)]-\ (23) 



and a structure factor dependent distinct hydrodynamic 
function part, 

sm{y' 



H%)\s,-^l dy 



y' 



[i + <f>s^„{y')r' 

/'^rfM(l-M')[5(|q-q'|)-l], (24) 



where /i is the cosine of the angle between q and q' [42] . 
The function S^g (y) , which should not be confused with 
the static structure factor, is given in [15, 42] as an in- 
finite sum of wavenumber-dependent contributions with 
inter-related scalar coefficients 7q"\ n = ... 00. Numer- 

(n) 

ical results for 7q obtained from a computation trun- 
cated at 77 = 5 have been given in Table 1 of the original 
paper by Beenakker and Mazur [15]. Taking advantage of 
the nowadays available computing power, we have been 
able to extend these earlier computations to more terms 
with truncations at ri = 10 and 15. However, our more 
accurate results for 7q"^ differ from the original results 
by Beenakker and Mazur by no more than 3%, and the 
differences in H'^{y)\s-f, ds{4>)\sf, and iioolsj are negligible 
for all practical purposes. 

The high-frequency limiting viscosity in the second- 
order 57-scheme is given by [17] 



Voo 

m 

Xo- 
A2 = 



1 



Ao + A2 ' 



(25) 



1 + 0-^70 



(2) 



1-^0 +§^0^ + 0(0^)^26) 



300 
An 



Ao7o 



^ jf(y/2) [S{y/<T)-l] 
1 + 057^ (77/2) ■ 



(27) 



where 7^^) ^ ^(^) ^ i + 167/840 -^0(0^). 
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Insertion of the low-concentration form, S^^{y) = 
1 — 240ji(y)/2/ + 0(^^), of the static structure factor of 
neutral hard spheres into Eqs. (23), (24), and (25) gives, 
after a straightforward calculation, the first-order virial 
expansion results 







do 






S-y 








S-y 



+ o(r) = 1 

+ 0(0') - 1 

+ 0(</.3) = 1 + ^0 + 

« 1 



131 
"56" 
411 
"56" 
5 
2 
5 
2 



(f>Kl- 2.339</>, (28) 

w 1 - 7.339(/), (29) 
1255 ,o 



168 



7.47?i^ 



(30) 



predicted by the scheme. The magnitudes of the lin- 
ear virial cocfEcients of d^^/do and K^^, and of the 
quadratic coefficient of tj}^ /rjQ, overestimate the precise 
values given in Eqs. (20), (21), and (22) by 28%, 12%, 
and 5% respectively. The effect of His on these quan- 
tities on the pair-level is thus overestimated by the ^7 
scheme. Clearly, the PA scheme is the method of choice 
when very dilute systems are considered. We note fur- 
ther that Hf^{q ^ 0) = -5(j) + 0{(j)'^) for hard spheres, 
a result quite close to the exact result of H'^{q — > 0) = 
-4.7140-1- 0(0'). This indicates that the zeroth-order 
(57 scheme is in general a better approximation for the 
distinct part, H'^{q), of the hydrodynamic function than 
for its self-part ds- 

Interestingly enough, the first-order in result for 
H''^{q -> 0) for hard-spheres predicted by the (57-scheme, 
is identical to the one obtained from the Rotne-Prager 
(RP) approximation of the His, where only the leading 
order monopole and dipole terms in the 1/x expansion of 
fifj are retained. For hard spheres, the first-order virial 
expansion result, H[q,n) = 1 — 1.350, for the principal 
peak height of H{q) remains valid to high accuracy up to 
the volume fraction 0/ = 0.494 at freezing []_4], whereas 
in the RP approximation, peak values of Hiji) larger than 
one are predicted. The main reason for this failure of the 
RP approximation lies in its prediction of ds = do at 
all 0, whereas the actual ds of hard spheres is signifi- 
cantly decreasing with increasing 0, down to the value 
(iP(0/) « 0.2do at freezing. 

In low-salinity charge-stabilized systems at low concen- 
trations, where qm<y w 27ran'^/^ <C 1, the (57-schemc gives 
predictions for H'^{q) close to those by the RP approx- 
imation. Indeed, these are precisely the systems where 
the RP approximation can be expected to perform well, 
explaining in part the overall success of the scheme in 
making reliable predictions for the distinct part of H{q) 
of charge-stabilized systems. 



C. Self- part corrected (57-scheme 

The key observation regarding the zeroth-order 6"f ex- 
pression for H'^(q), which depends on S{q) only, is that 



it gives overall good results both for neutral and charged 
Yukawa-type spheres. In contrast, the zeroth-order Sj 
expansion for ds in Eq. (23) depends on only, indepen- 
dent of the employed pair-potential. Comparison with 
ASD simulation results [14], and experimental data for 
ds for charged colloids [!), 4.']], show that Eq. (23) is a de- 
cent approximation of ds for neutral hard spheres only. 

The self-diffusion coefRcient, ds{(j>), of charged spheres 
is in fact larger than the one for neutral spheres at the 
same [44], since for the latter near contact configura- 
tions are more likely. Using leading-order far-field mobil- 
ities applicable to strongly charged colloids characterized 
by qm oc (j)^^^, one finds for < 0.1 a power- law depen- 
dence of ds according to ds/do ~ 1 — atcj)^/'^ [2(S, 44-47], 
differing qualitatively from the regular hard-sphere virial 
result in Eq. (20). The coefficient at ~ 2.5 — 2.9 in the 
fractional power law varies to a certain extent with the 
particle size and charge. 

For suspensions of strongly charged spheres, where 
< 0.15, it has been shown [9] that the PA result for ds 
in Eq. (14) is in better agreement with ASD simulation 
results, and experimental data, than the corresponding 
(57-scheme result based on Eq. (23). For larger concentra- 
tions > 0.15, the PA scheme overestimates the slowing 
hydrodynamic influence on ds, since it does not account 
for the shielding of the His in pairs of particles by other 
intervening particles. Hydrodynamic shielding is a many- 
body effect which lowers the strength but not the range of 
the His. It should be distinguished from the screening of 
His by spatially fixed obstacles or boundaries which ab- 
sorb momentum from the fluid, thereby causing a faster 
than 1/r decay of the flow perturbation created by a 
point-like force. The neglect of hydrodynamic shielding 
(i.e., three-body and higher-order His) by the PA scheme 
is more consequential for the sedimentation coefficient K 
than for ds- To the former, the PA scheme is applicable 
to decent accuracy only up to w 0.1 [9], whereas for 
larger 0, the coefficient K becomes increasingly under- 
estimated. The coefficient ds, on the other hand, is less 
sensitive to the neglect of higher-order HI contributions 
than K or dc, since the leading-order far- held contribu- 
tions to Xii{x) and yii{x) are of 0(x~^), i.e., of shorter 
range than the leading-order 0(x~^) contributions to K . 

As a simple improvement over the zeroth-order 5^ 
scheme for the H{q) of charged particles, which preserves 
its analytic simplicity, we therefore use the self-part cor- 
rected expression 



H{y)\ 



Sic. 



ds_ 
do 



PA 



57' 



(31) 



with ds according to Eq. (14) and H'^{q) according to 
Eq. (24), bearing in mind that [ds/do]pA becomes less 
reliable for > 0.15. 

For the limiting case of neutral hard spheres at larger 
0, it is therefore preferential to use in place of [ds/do]pA 
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the accurate expression, 

HS 



do 



1- 1.83150(1 + 0.11950 -0.70(/)2) (32) 



which, for < 0.5, agrees well with ASD [14] and hydro- 
dynamic force multipole [48] results, with an accuracy 
better than 3%. Note that the linear and quadratic or- 
der coefficients in Eq. (32) have been selected identical 
to the numerically precise values —1.8315 and —0.219 = 
— 1.8315 X 0.1195, for the respective virial coefficients 
given in Ref. [39]. We have determined the cubic co- 
efficient in Eq. (32) from a best fit to recent simulation 
results in [48] and [14]. The coefficient 0.70 differs some- 
what from the coefficient 0.65 in [48], where simulation 
results only up to < 0.45 were considered. 

A self-part corrected version of the scheme for H{q) 
was used already in earlier applications, where ds was 
considered simply as an adjustable parameter [42], or in 
more recent work determined using elaborate ASD sim- 
ulations [14]. For practical purposes, however, it is far 
more convenient to use the analytic dg corrections in 
Eqs. (14), (31) and (32). In the present work, numer- 
ous ASD simulation results for H{q), ds, and rjoo have 
been generated to provide precise benchmarks for assess- 
ing the accuracy of the proposed self-part corrected 5^ 
scheme. 

While the 5^ scheme for H{q) has been used already 
for charge-stabilized colloids, to our knowledge the ap- 
plication of the ^7 scheme for rj^o in Eqs. (25)-(27) was 
restricted so far to colloidal hard spheres, where the pre- 
dicted values for ?7oo(</') are in good agreement, for < 
0.4, with experiments and simulation data. In the fol- 
lowing section, we are going to assess the performance of 
the (57-scheme expression for the rjoo of charged particles 
by comparison with ASD simulation data. Our analysis 
shows that the second-order contribution, A2(0) < 0, 
to iioo/vo is only weakly dependent on the shape of the 
static structure factor. Moreover, the zeroth-order con- 
tribution, Ao(0) > 0, which is only dependent on 0, dom- 
inates for small 4> the contribution A2 in magnitude. As 
a consequence, the (57-predicted values for 7700 change 
only slightly when going from neutral to charged parti- 
cles, whereas simulations, and experiments [49], reveal 
significantly smaller viscosity values in particular at low 
salinities. Thus, the 57-schcmc result in Eqs. (25)-(27) 
applies to neutral hard spheres only. However, for the in- 
teresting case of low-salinity systems, where the viscosity 
differences to neutral spheres at equal concentration are 
largest, the 5^ scheme can be modified (corrected) in an 
ad-hoc way, according to 



(57c, 



1 



1 



Aq Ao + A2 



(33) 



The motivation for this correction follows from the PA ex- 
pression in Eq. (18): For a low-salinity system of strongly 
repelling particles, one has the scaling oc for 
the peak position, r™, of the rdf, and g{r < r,„) « 0. 



Since J{x) is of 0{x~^), for these systems the integral in 
Eq. (18) is of 0(0'^). Hence, to quadratic order in cj), rjoo is 
determined basically by the microstructure-independent 
contribution, 1 -I- 2.50(1 -|- 0), to Eq. (18). 

In Eq. (33) , we correct approximately for this limiting 
behavior of 7700 by subtracting the structure-independent 
"self part", I/Aq, from [ijoo/vols^^ which renders the 
remainder of 0(0^) small, while adding the term 1 -|- 
2.50(1 + 0). As we will show in the following, the so- 
corrected (^7 scheme is in very good agreement with the 
ASD viscosity data of low-salinity systems, even up to 
the freezing transition concentration. We point out here 
that Eq. (33) is restricted in its applicability to the low- 
salinity regime of strongly repelling particles. In con- 
trast, the dj-corrected 57-scheme for H{q) in Eq. (31) 
applies to HSY systems for any set of system parame- 
ters {7, fc, 0}, provided < 0.15, including the crossover 
regime from neutral to dcionizcd, highly charged particle 
systems. For neutral hard spheres, the (uncorrected) ^7 
scheme for ?yoo performs quite well. 

The design of a simple, corrected ^7 scheme which op- 
erates well for arbitrary {7, fc, 0}, including systems of 
intermediate salinity, is obstructed by the limited sep- 
arability of Aq and A2, and by significant many-body 
HI contributions for more concentrated systems of nearly 
hard-sphere-like particles at high salinity. For small val- 
ues < 0.1, the PA method can be used to produce 
reliable predictions of ryoo; for arbitrary salinities. 



D. Accelerated Stokesian Dynamics simulations 

The simulation data for the H{q) and rjoo in HSY 
systems explored in this work, have been generated us- 
ing an accelerated Stokesian Dynamics (ASD) simulation 
code. The details of the simulation method have been 
explained in Ref. [32]. It allows to simulate short-time 
properties of a larger number of spheres, typically up to 
N = 1000, placed in a periodically replicated simulation 
box, allowing for improved statistics. Since short-time 
properties are obtained from single-time equilibrium av- 
erages, we have used equilibrium configurations gener- 
ated using a Monte-Carlo simulation method for charged 
spheres and a Molecular Dynamics algorithm for neu- 
tral hard spheres, with the many-sphere His accounted 
for using the ASD scheme. The computed hydrodynamic 
function, HN^q), shows a strong system-size dependence, 
even when N is not small. We therefore extrapolate 
Hjsi^q) to the thermodynamic limit using the finite-size 
scaling correction [50, 51], 



H{q) = HN{q) + 1.76Siq) 



Vo 



^700(0) 



(0/iV) 



1/3 



(34) 



which, for g — >■ and q 00, includes the finite-size cor- 
rections for K and ds, respectively. This finite-size cor- 
rection formula was initially proposed by Ladd for hard 
spheres [50, 52, 53], and has been subsequently applied 
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also to charged spheres [14] and solvent-pcrmcablc parti- 
cles [3-3]. As pointed out by Ladd [50], and explained by 
Mo and Sangani [51], 7700 is not critically dependent on N 
so that finite size scaling extrapolation to 00 is not 
needed (see also [.'>4]). The simulation results discussed 
in following arc obtained from averaging over 2000 con- 
figurations, for systems of typically N ~ 512 particles. 



V. RESULTS 
A. Diffusion properties of charged particles 

In Fig. 1, theoretical and simulation results for H{q) 
are shown for nine different systems with common system 
parameters Lb = 5.617 nm, a = 200 nm, and Z = 100, 
representative of suspensions of highly charged colloidal 
spheres in an organic solvent. The parameters (j) and 
are varied, and assume values consisting of all permuta- 
tions of 0.055, 0.105, and 0.15 and = lO"'', 10-^ 
and 10~^M. To facilitate the comparison of the different 
systems, identical axes scales are used in all nine panels, 
ordered with respect to 0, which increases from top to 
bottom, and Tig, decreasing from left to right. Thus, the 
strength of the interparticlc correlations increases from 
left to right, and from top to bottom. Fig. 1 serves for 
analyzing the accuracy of the corrected and uncorrected 
6^ schemes, and of the PA scheme, in comparison to our 
finite-size corrected ASD results for H{q). 

The results for H{q) obtained by all methods described 
in Sec. IV, are included in Fig. 1. Open symbols repre- 
sent ASD simulation data, black dashed curves our PA 
results, blue dotted curves the zeroth-order ^7 scheme, 
and red solid lines the self-part corrected Sj scheme pre- 
dictions. The black and orange horizontal lines in each 
panel mark the reduced short-time self-diffusion coeffi- 
cient, dg/do = H{q — s> 00), obtained from the ASD sim- 
ulations and the PA scheme, respectively. As the only 
input to the three analytic schemes, S{q) and g{r) of 
each system are obtained using our MPB-RMSA code. 
Except for (/> = 0.105, the MPB-RMSA resuhs for S{q) 
and g{r) have been shown already in Figs. 3 and 4 of [18], 
and are therefore not included here. 

The rightmost column of panels in Fig. 1 presents re- 
sults for three systems of strongly charged particles with 
a very low residual, but experimentally still accessible, 
salt content. In the most concentrated system in panel 
(i), where (f> ~ 0.15 and = 10^® M, a structure factor 
peak value S{q„i) « 2.8 is attained according to the MC 
simulations. The very same peak value is predicted by 
the MPB-RMSA and RY integral equation schemes. Ac- 
cording to the empirical Hansen- Verlet freezing rule, the 
system in panel (i) is pretty close to the freezing transi- 
tion point [54-56]. 

The screening parameter k defined in Eq. (3) assumes 
rather low values of 2.67, 3.24, and 3.68 for the systems 
in panels (c), (f), and (i), with k1/k1 = 1.1, 2.1. and 
3.0, respectively. The relatively large values for k^/kl 



in these systems indicate that salt microions contribute 
little to the electrostatic screening, which instead is dom- 
inated by colloid-surface released counterions. Different 
from neutral hard spheres, where H{qm) decreases lin- 
early with increasing </>, the hydrodynamic function peak 
heights of the three low-salinity systems depends non- 
monotonically on 0, with ASD values H{qm) = 1.13, 
1.17, and 1.15 for the systems in panels (c), (f), and 
(i), respectively. Such a non-monotonic </)-dependencc of 
H{qm) is typical for low-salinity systems, as discussed 
in [14] and [57]. The ASD results for the reduced self- 
diffusion coefficient, ds/do, and the corrected (57-scheme 
results for the sedimentation coefficient K in panels (c), 

(f) , and (i), follow closely the concentration-scaling pre- 
dictions dg/do = 1 — atcj)^^^ and K = 1 — agcj)^/^ given in 
[!)], with at = 2.63 and with Os = 1.44. The fractional 
exponents 4/3 and 1/3 can be traced back to the (j?-/^ 
scaling in low-salinity systems, of the structure factor 
peak position, q„i. Note here that H{q) has its principal 
peak at a wavenumber practically identical to the peak 
position of S{q). 

We proceed in our discussion with the systems in 
the leftmost column of panels. The salt concentration. 
Us = lO"'* M, of these systems is so large that the neutral 
hard-sphere (HS) limit is practically reached. The com- 
parison with the H{q) of genuine neutral hard spheres at 
volume fractions equal to those in panels (a), (d), and 

(g) ; shows relative differences of less than 6% in all three 
cases. The proximity to genuine neutral hard-sphere sys- 
tems is manifest also in the large values, k = 18.5, 18.6, 
and 18.7, of the screening parameter, and in the small 
ratios k^/k^ = 0.01,0.02, and 0.03, for the systems in 
panels (a), (d), and (g), respectively. 

For the smallest considered concentration (f> = 0.055 
(top row of panels in Fig. 1), the differences in the re- 
spective H{q) predicted by the analytic methods and the 
ASD simulations are very small. Since the PA scheme be- 
comes exact at low 0, this illustrates that, despite their 
overall inaccuracies, the self-part corrected, and even the 
uncorrected (57-scheme, can be used to obtain good esti- 
mates of H(q) also for more dilute suspensions. 

With increasing pronounced differences are observed 
in Fig. 1 between the PA-scheme and ASD results for 
H{q). This reflects the expected failure of the PA 
scheme in concentrated suspensions, where three-body 
and higher-order HI contributions become influential. 
The deviations of the PA-scheme results from the pre- 
cise simulation data are most pronounced for the peak 
value H{qm), which is overestimated by the PA scheme, 
and in the sedimentation coefficient, K = H{q — 0), 
which is underestimated. In fact, for the system in panel 
(i), the PA prediction for K is just barely larger than 
zero, turning to unphysical negative values when the vol- 
ume fraction surpasses = 0.154 at a fixed Ug = 10^^ 
M. However, the PA-scheme values for dg/do remain in 
very good agreement with the ASD results, with a rela- 
tive deviation of less than 3.5% even at = 0.15. The 
values for dg/do predicted by the uncorrected Sj scheme 
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Fig. 1. Hydrodynamic function, H{q), for charge-stabilized suspensions of volume fractions, (j), and salt concentrations, n^, 
as indicated in each panel. The panels are ordered with respect to (f>, which increases from top to bottom, and ria, which 
decreases from left to right. Open symbols: ASD simulation data. Black dashed, blue dotted, and red solid curves: PA-scheme, 
(57-scheme, and self-part corrected 57-scheme results, respectively. Horizontal black and orange lines mark the values for ds/do 
obtained from ASD simulation and PA-scheme calculations, respectively. Common system parameters are Lb = 5.617 nm, 
a = 200 nm, and Z = 100. 



are generally in less good agreement with the simulation 
data, clearly revealed in Fig. 1 by the large-g offset of the 
corresponding H{q). 

The self-part corrected (57-scheme results for in 
Fig. 1 (red solid lines) illustrate that this hybrid scheme 
combines the good accuracy of the PA scheme regarding 
ds, and of the 5^ scheme regarding H'^{q). Indeed, the 
corrected ^7-scheme results for H(q) are in overall good 



agreement with the ASD simulation data for all consid- 
ered systems, with the largest deviation of 6% for H{qm) 
observed in panel (i). 

In closing our discussion of Fig. 1, a short comment is 
in order regarding the computational cost caused by the 
considered methods of computing H{q). The fast and 
accurate evaluation of S{q) and g{r) by the MPB-RMSA 
method, in combination with the easily evaluable inte- 



10 




Fig. 2. Reduced short-time sedimentation coefHcient, K , 
of neutral colloidal hard spheres. Open circles: Hydrody- 
namic force multipole simulation data by Abade et al. [33]. 
Open Squares: Lattice-Boltzmann simulation data by Segre 
et al. [o2]. Black dashed line: PA-scheme result. Dashed- 
dotted red line: uncorrected Sy-scheme result. Dashed orange 
line: self-part corrected 57-scheme result, with ds/dg taken 
from the PA-scheme. Solid black line: self-part corrected S"/- 
scheme result, with dg/dg according to Eq. (32). Solid blue 
line: second-order virial result K^^ = 1 - 6.546(?!) -I- 21.9180^ 
[40]. The static structure factor input was calculated using 
the analytic Percus-Yevick solution. 

grals in Eqs. (14), (15), and (24), has allowed us to imple- 
ment a convenient graphical user interface code, running 
on a standard desktop PC. Using this code, MPB-RMSA 
results for S{q) and 5(r), and PA-, 5^-, and self-part cor- 
rected (57-scheme results for H{q), are obtained in less 
than 1 second of cpu time, for a given set of input pa- 
rameters {Lb, (T, .Z, ris, 0}. Thus, all curves depicted in 
Fig. 1, except for the ASD simulation data, have been 
obtained altogether in less than a minute on a standard 
desktop PC. In comparison, the computation of just one 
of the computer simulation curves in Fig. 1 required on a 
standard desktop PC typically 5 hours of cpu time for 
generating 2000 equilibrated configurations (using our 
MC method) and approximately 8 hours of cpu time for 
computing H{q) with the ASD scheme. The overall ac- 
curacy and fast performance of the hybrid ^7 scheme in 
Eqs. (31)-(33) make this scheme well-suited for the real- 
time fitting, even of large sets of experimentally recorded 
data for H{q) and D{q) ["iS]. 

B. Hybrid ^7 scheme applied to neutral hard 
spheres 

The main virtue of the c?s-corrected (zeroth-order) ^7 
scheme lies in its good applicability to charge-stabilized 
systems. However, it is interesting to assess in more de- 
tail its performance in the limiting case of neutral hard 
spheres, in particular when the values of H(q) at q = 
and q = qm arc considered. Recall for hard spheres that 



the accurate expression for df^ in Eq. (32) should be 
preferentially used instead of the approximate PA result 
for larger (p > 0.15. For neutral spheres, higher-order HI 
contributions to ds begin to matter at somewhat smaller 
concentrations than for charge-stabilized particles, where 
near-contact configurations are unlikely. 

In Fig. 2, numerically precise Lattice Boltzmann [52] 
and hydrodynamic force multipole simulation results [33] 
for K^^{(j)) are compared with the predictions of all con- 
sidered analytical schemes. For cj) < 0.35, the uncor- 
rected 6j scheme underestimates the simulation data, 
showing the opposite trend of a slight overestimation for 
(f) 0.35. The corrected 5^ scheme with dg-input accord- 
ing to Eq. (32), on the other hand, is in excellent agree- 
ment with the simulation data up to ~ 0.4, reflecting 
the accuracy of the 57-scheme predictions for H'^[q) also 
for neutral spheres. 

For large volume fractions (f> > 0.4, however, the dis- 
tinct part, K^^ — df^/do, of the sedimentation coeffi- 
cient is considerably underestimated by the ^7 scheme, 
to such an extent that the self-part corrected Sj scheme 
prediction for K^^ assumes unphysical negative values 
for (j) > 0.45. Up to w 0.2, the corrected 57-schcme pre- 
diction for K^^, with df^ obtained by the PA scheme, 
lies closer to the simulation data than the uncorrected 
(57-scheme result. This can be explained by the pre- 
cise account of (two-body) lubrication effects in the PA- 
scheme, which are not included in the uncorrected ^7 
scheme. At larger concentrations, however, the corrected 
6 J scheme, with PA input for ds, increasingly underesti- 
mates the sedimentation coefficient up to the point that, 
for > 0.31, unphysically negative values for K^^ are at- 
tained. This is a consequence of the already noted many- 
sphere hydrodynamic shielding effect, disregarded in the 
PA scheme, which lowers the strength of the His with- 
out reducing their range, leading to a larger self-diffusion 
coefficient than predicted on basis of hydrodynamic pair- 
interactions alone. The neglect of shielding effects by the 
PA scheme, both in the self- and distinct parts of K^^, 
is the reason for the crossover of the PA curve of K^^ to 
negative values already at (/> « 0.21. 

Regarding again Fig. 2, we finally note that the second- 
order virial resuh, K^^ = 1 - 6. 546(/)-h 21.91802 + 0(0^)^ 
derived in [4(J] ceases to be applicable for (j> > 0.15, where 
its curve bends up to larger values. This is the reason why 
this second-order virial result cannot be used, different 
from the corresponding virial results for df^, f]^ , and 
H^^{qm), to construct analytic extrapolation formulas, 
valid for all concentrations up to the freezing transition. 

We proceed by discussing the concentration depen- 
dence of the peak value, H^^{q,n), of the hydrodynamic 
function of neutral spheres. As noted in Sec. Ill, H(qm) 
is related to the short-time cage diffusion coefficient, 
dcgc = doH{qm)/ S{qm), characterizing the initial decay 
rate of density fluctuations of wavelength equal to the 
next-neighbor cage size. Fig. 3 displays the decline of 
the hard-sphere H(qm) with increasing (j). To excellent 
accuracy up to the freezing volume fraction, this decline 
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Fig. 3. Hydrodynamic function peak value, H[qm), of neutral 
hard spheres. Open circles: Hydrodynamic force multipole 
simulation results by Abade et al. [.>.!]. Open diamonds: ASD 
simulation data [11]. Black dashed line and dashed-dotted 
red line: PA-scheme and uncorrected J7-scheme results, re- 
spectively. Dashed orange line: self-part corrected 57-scheme 
results, with da/do taken from PA-scheme calculations. Solid 
black line: self-part corrected 5^-scheme results, with ds/do 
according to Eq. (32). Dotted curve in magenta: 1 — 1.350. 
The static structure factor input was obtained using the an- 
alytic Percus-Yevick solution. 



is described by the first-order virial result [14] 

i7HS(g„j = 1 _ 1.350. (35) 

Indeed, all the depicted ASD [14] and hydrodynamic 
force multipole [33] values for H{qm) follow this line, in- 
dicating that, for a so far unknown reason, all higher or- 
der virial contributions cancel out. According to Fig. 3, 
the uncorrected scheme significantly underestimates 
H^^{qm) for 4> ^ 0.35, overestimating it instead for 
(j) > 0.4. In contrast, the corrected ^7-schemc result 
with the precise d^^ according to Eq. (32), is distinctly 
more accurate in that it only very slightly underestimates 
the linear decay in Eq. (35) for (p < 0.4. Moreover, for 
(j) > 0.4, the positive- valued deviations from 1 — 1.35(f> 
are substantially smaller than those of the uncorrected 
(57 scheme. 

The corrected J7-scheme prediction for H^^{qm), with 
ds calculated using the PA scheme, is a decent approxi- 
mation up to (/) < 0.15. Its bending over to smaller val- 
ues occurs for H^^{qm) at somewhat larger cj) than in the 
sedimentation case, indicating that (57-scheme results for 
H'^{q) arc more accurate aX q = q^ than at q w 0. The 
curve for H^^{qm) predicted by the PA scheme bends 
over to larger values at a concentration (j) « 0.37 vastly 
beyond its range (0 < 0.1) of applicability. The total 
neglect in the PA scheme of many-body His beyond the 
pair level implies, at larger 0, an underestimation of dg, 
but to a larger extent an overestimation of H'^{qm)- As a 
net result, H^^{qm) at large (j) is strongly overestimated 
by the PA scheme. 
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Fig. 4. Reduced high-frequency viscosity, rioo/rjo, as a func- 
tion of (j), for a suspension of neutral hard spheres (HS, in 
black), and two deionized charged-sphere suspensions (CS-1 
and CS-2, in red). The leading-order Einstein contribution, 
1 + 2.50, is subtracted off to expose the differences. Symbols: 
ASD simulation results. Dashed, dotted and solid lines: PA- 
, S-y-, and self-part corrected 57-scheme results, respectively. 
All analytic schemes use the MPB-RMSA S{q) as input. The 
CS-1 viscosity results represented by red filled circles are ASD 
data for Ls = 5.617 nm, a = 200 nm, and Z = 100. The ASD 
data for the more weakly charged, smaller particles of system 
CS-2, with Lb ~ 0.71 nm, cr = 50 nm, and Z — 70 are in- 
dicated by red diamonds filled in blue. The parameters of 
system CS-2 have been used in the analytic calculations. The 
inset magnifies the details at lower 0. 

In summarizing our discussion of hard-sphere systems, 
the key message conveyed by Figs. 2 and 3 is that the 
corrected ^7 scheme, with ds according to Eq. (32), de- 
scribes H^^{q) quite precisely for < 0.4. It can be 
applied to reasonable accuracy even to larger values, 
with the exception of small q values. 



C. High-frequency viscosity 

In the following, we compare viscosity results obtained 
by the various methods described in Sec. IV, for the two 
limiting cases of deionized (low-salinity) charged-sphere 
and neutral hard-sphere suspensions. Results for systems 
with intermediate added salt are bracketed by these two 
limiting cases. 

In Fig. 4, we display our viscosity results for neutral 
hard spheres (HS) with those for two deionized suspen- 
sions of highly charged spheres (CS) where Ug = 0. Re- 
sults by all the methods in Sec. IV are shown. We point 
out that, in addition to the CS system of Fig. 1 with pa- 
rameters Lb = 5.617 nm, a = 200 nm, and Z = 100 (re- 
ferred to here as system CS-1), whose ASD results for t]^ 
in Fig. 4 are indicated by filled red circles, we additionally 
show viscosity results for another zero-salt system where 
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Lb = 0.71 nm, cr = 50 nm, and Z = 70 (labeled CS-2), 
whose ASD data for rjoo are represented by red diamonds 
filled in blue. The reason for including in Fig. 4 results 
for two different deionized systems, is that system CS- 

1 freezes at « 0.15, whereas systems CS-2 stays fluid 
up to ^ 0.3, allowing us to test the predictions of our 
analytical methods in a more extended volume fraction 
range. The ASD simulation data for ?7oo(0) for systems 
CS-1 and CS-2 merge continuously, overlapping nearly 
perfectly within 0.1 < (j> < 0.15. This indicates that 
the limiting behavior of rjoa for highly correlated charged 
spheres is practically reached in both systems. There- 
fore, the depicted CS-results for rjoo were calculated in 
the analytic schemes using the parameters of system CS- 

2 only. The deviations in rjoo for the two CS systems are 
minuscule in all considered analytic schemes. 

Consider first the performance of the PA scheme. By 
its definition, the PA result for 7700 is in very good agree- 
ment with the ASD simulation data at low 0, but the 
agreement becomes poorer with increasing volume frac- 
tion. While the deviations between the ASD and PA 
results for low-salt systems are very small up to (f> < 0.2, 
for neutral spheres significant differences are visible al- 
ready for > 0.1. These differences originate from 
the fact that in charged- sphere suspensions, near-contact 
configurations are disfavored by the electric repulsion, 
i.e., charged-sphere systems are hydrodynamically more 
dilute than neutral sphere suspensions. Since higher- 
order HI effects on 7700 in low-salinity charged systems 
are weaker, for many such systems (including system CS- 
1), which freeze already at < 0.2, the accuracy of the 
PA-scheme for t]oo is sufficiently good in the whole fluid 
regime. Regarding H(q), however, the PA-scheme pre- 
dictions for charged spheres deviate significantly from the 
ASD data already at = 0.105 (see again Fig. 1). At 
larger cj), and in contrast to the ASD data, the PA scheme 
predicts only a slight enlargement of 77^0 in going from 
charged to neutral spheres. The distinctly larger values 
of 7700 for concentrated hard-sphere suspensions are thus 
mostly due to near-field, many-body His which enlarge 
the viscous dissipation. Overall, however, 7700 is rather 
insensitive to the range of the pair potential, at least in 
comparison to the static (zero frequency) viscosity which 
for concentrated systems can become very much larger 
than 7700 [59]. 

The self-part modified Sj scheme for 7/00, defined by 
Eq. (33), agrees overall very well with the ASD data 
for charged spheres in the whole fluid-state concentra- 
tion regime. Small deviations from the simulation data 
are noticed at low (f> values only. Regarding neutral hard 
spheres, a similar observation applies to the unmodified 
second-order S"f scheme, which describes the ASD simu- 
lation data quite well up to « 0.4. The slight overesti- 
mation of 77^^ at lower cj) can be attributed to the non- 
exact treatment of two-body HI contributions by the ^7 
scheme. 

Overall, the high-frequency viscosity of charged-sphere 
systems at low salinity is well captured by the modified 



Sj-scheme in Eq. (33), and for neutral hard spheres by 
the unmodified Sj scheme (up to (/) « 0.4). Different 
from the self-part corrected ^7 scheme for H{q), which 
makes reliable predictions for arbitrary salinities, the 
modified scheme for rjoo applies to low-salinity sys- 
tems only, and the unmodified scheme only to neutral 
hard spheres. The reasons for this have been discussed 
already in subsection IV C. 

Before closing our discussion of rjoo , it is of interest to 
compare the numerical efforts required by the employed 
methods. The computation of the 45 ASD data points 
for neutral and charged spheres included in Fig. 4 re- 
quired about 500 hours of cpu time on a modern desktop 
PC. This large time investment should be compared to 
the few minutes computation time on a comparable PC 
which were required for the results by all considered an- 
alytic schemes, amounting to more than one thousand 
data points on a dense mesh of values. 

VI. RELATION BETWEEN VISCOSITY AND 
DIFFUSION PROPERTIES 

Having quite accurate analytic methods for short-time 
properties at our disposal, we are in the position to 
analyze possible relations between these properties in 
the whole fluid-state concentration regime. Specifically, 
we want to test the validity of two generalized Stokes- 
Einstcin (GSE) relations, 

^x^.l, (36) 
ao Va 

for D*[4>) = ds{(i>) and D*{(f>) = dcgc(0)- In addition, 
we probe the validity of the Kholodenko-Douglas GSE 
(KD-GSE) relation [60], 

^x!^xV^(,^O,0)«l, (37) 

between the collective diffusion coefficient, dc = 
K/S{q — > 0), 7^00, and the square root of the isothermal 
osmotic compressibility given by S{q — ?> 0). In particular 
the KD-GSE relation has been used in various biophysi- 
cal and soft matter studies [Gl-64]. 

All three considered GSE relations are exact at = 
only. The approximate validity of a GSE relation in con- 
centrated systems is an important issue in microrheolog- 
ical studies, since this allows to infer a rheological prop- 
erty more easily from a diffusion measurement. For test- 
ing the GSE relations in Eqs. (36) and (37), we consider 
here again the two limiting HSY cases of a low-salinity 
charge-stabilized system and neutral spheres, since the 
differences in the respective short-time dynamic proper- 
ties are here largest. 

For a precise test of the GSE relations in the case of 
hard spheres, we take advantage of simple analytic ex- 
pressions available for all short-time properties appear- 
ing in Eqs. (36) and Eq. (37), with the exception of K, 
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Fig. 5. Test of the two GSE relations in Eq. (36), relating rjoo 
to ds and dcge, respectively. Solid lines: (ds/do) x (rjoo/vo) vs. 
(j). Dashed lines; (dcgc/do) x {rjoo/ilo) vs. 0. Black curves are 
hard-sphere results obtained from Eqs. (32), (35), (38), and 
(39). Red and blue curves are corrected (57-scheme results, 
using as input the MPB-RMSA results for S{q) of the low- 
salinity systems CS-1 and CS-2, for — 10~^ M. 



for which we use the quite accurate self-part corrected 
i57-scheme result depicted in Fig. 2. The analytic expres- 
sions for hard spheres, which apply to excellent accuracy 
up to = 0.5, are Eq. (35) for H^^{qm), Eq. (32) for 
df^, and the generalized Saito-type expression for rj^ 
[34], 



n, 



HS 



1 



1 + 5* 



2^1-0(1 + 5)' 



(38) 



where 5= l.OOl(/)+O.9502-2.15(/)3. Moreover, we use the 
precise formula [5!)] for the structure factor peak height 



:<HS 



(g™)«l + O.6440.g^^(x-l+), 



(39) 



where g^'^ix = 1+) = (1 - O.50)/(l-0)^ is the 
Carnahan-Starling contact value for the hard-sphere rdf. 
For S^^{q — > 0) in Eq. (37), we employ the Carnahan- 
Starling equation of state [22] . 

In testing the GSE relations in Eqs. (36) and (37) for 
low-salinity systems, for the diffusion properties we use 
the self-part corrected 6j scheme, with ds calculated by 
the PA scheme. For 7700 we use the corrected (S7 scheme 
according to Eq. (33). 

The validity of the two GSE relations in Eq. (36) is 
examined in Fig. 5. A valid GSE relation is reflected by 
a horizontal line of unit height. For neutral hard spheres 
(black lines), the product (d^^/do) x /ijq) is well 
approximated, for all displayed 0, by its first-order in (jj 
expansion given by 1 4- 0.670, showing a more than 20% 
violation of this GSE relation for d^^ when cfi > 0.3. 

Different from d^^ , the GSE scaling for d^^^ is approx- 
imately satisfied with a maximal deviation from one of 
8%. Thus, for hard-sphere like colloidal particles avail- 
able only in amounts too small for a mechanical rhco- 
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Fig. 6. Test of Kholodenko-Douglas GSE relation in Eq. (37). 
Black curves: neutral hard-sphere results based on precise an- 
alytic expressions for 77^^, S'^^(q — > 0), and K^^ calculated 
in self-part corrected 57-scheme, with ds taken from Eq. (32). 
Red (dashed) and blue (dashed-dotted) curves are self-part 
corrected (57-scheme results for the low-salinity systems CS-1 
and CS-2 (using — 10~^ M), respectively, with ds calcu- 
lated using the PA scheme, and 5(g) using the MPB-RMSA 
method. 



logical experiment, one can determine r]rx, approximately 
from a dynamic scattering experiment measuring D(qm)- 

According to Fig. 5, in low-salinity systems of charged 
particles, the GSE-relation for ds is overall of similar 
accuracy as that for hard spheres, although the devia- 
tions from one are larger at smaller (j). The curves for 
(ds/do) X (7700 /?7o) obtained for the two low-salinity sys- 
tems coincide practically for (j) < 0.15. The downturn 
of these two curves at larger cfi > 0.18, indicated by the 
dotted curve continuations in Fig. 5, is not shared by the 
ASD simulation data (c.f.. Fig. 25 in [14]). This is an 
artifact of the PA scheme which, as discussed already in 
relation to Fig. 1, tends to underestimate ds at larger <j>. 

Different from neutral hard spheres, in low-salinity sys- 
tems the GSE relation for dcgc is manifestly violated 
already at very low (j)- The strong difference in the 
(dcgo/do) X ('7oo/'7o) curves for the two considered low- 
salinity systems is due to the different 0-dependence of 
their respective S{qm)- The pronounced decline of both 
curves at low (jj is mainly triggered by the sharp low-0 
rise of S{qm) in low-salinity systems. 

Kholodenko and Douglas [60] have proposed the GSE 
relation in Eq. (37) using mode-coupling theory like ar- 
guments. For neutral spheres at low 0, we can check 
this relation analytically using the numerically precise 
second-order virial expansion results for K^^ in [40], -q^ 
in [41], and S^^{q -> 0) given by the Carnahan-Starling 



equation of state. This leads to 



~1 — 



S'HS(q ^ 0,0) = 1-0.0460+1.371302+0(03), 

(40) 
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where the first-order virial coefficient is indeed close to 
zero. Using the analytic expressions for 77^^ in Eq. (38), 
and K^^ calculated in the self-part corrected (57-scheme 
with ds according to Eq. (32), we can test the KD-GSE 
relation in fact in the full fluid-state regime of neutral 
spheres. According to the inset in Fig. 6, the KD-GSE 
for neutral spheres is valid to decent accuracy up to < 
0.4, with a maximal positive-valued deviation from one 
at (/) w 0.3 of less than 18%. Strong negative-valued 
deviations from one are observed for 4> ^ 0.4, where the 
KD-GSE relation ceases to be applicable to neutral hard 
spheres. 

The corrected (57-scheme results included in Fig. 6 
demonstrate the striking violation of the KD-GSE re- 
lation, when low-salinity systems of charged particles are 
considered. A clear violation of this relation is observed 
for all concentrations cj) > 10^^ in the case of the low- 
salinity system CS-1, and for cf) > 10~^ in the case of 
system CS-2, where, for both systems, a residual salt 
concentration of Us = 10~^ M has been assumed. The 
maximal (positive- valued) violation of the KD-GSE rela- 
tion occurs roughly at a volume fraction where the dc{4') 
of charged spheres attains its maximum. The maximum 
in dc{4>), in turn, is the result of a competition, with 
increasing (f>, between decreasing compressibility and de- 
creasing sedimentation coefficient. The concentration at 
the peak of dc(0) is determined roughly from k'^{(f>) = k1 
[5]. The downturn of the dc'rioo\/S{q — > 0)/(do??o) curve 
at large 0, observable in Fig. 6 for neutral and charged 
spheres alike, is triggered by the large- decline both of 
K and r]oo- 

VII. CONCLUSIONS 

We have presented a comprehensive theoretical and 
computer simulation study of short-time dynamic prop- 
erties of colloidal spheres, interacting by a HSY pair po- 
tential. In this study, the accuracy of essentially analytic 
and very fast methods of calculating H{q), dg and rjoo has 
been explored in comparison with numerically expensive 
ASD simulation results, obtained as functions of volume 
fraction and added salt content. We have been particu- 
larly concerned with the low screening (low salinity) and 
infinite screening (i.e., neutral sphere) limits of the HSY 
model. 

The static pair functions S{q) and g{r), required as 
the only input to the PA and (self-part corrected) ^7 
schemes, have been determined using our recently de- 
veloped MPB-RMSA method. The PA scheme, which 
precisely (and only) accounts for two-body His, is for ar- 
bitrary salt concentration in excellent agreement with the 
simulation data for H{q), provided that < 0.1. With 
regard to the high-frequency viscosity, the PA-scheme 
predictions are in good agreement with the simulation 



data for volume fractions up to </> « 0.1 for neutral hard 
spheres, and up to (/> « 0.2 for strongly charged spheres 
in the weak screening (low salinity) regime. At larger cj), 
three-body and higher-order His become influential. 

The self-part corrected S'j scheme for H{q) is in good 
agreement with our ASD results for the hydrodynamic 
function of charged spheres at all considered </> and Hs 
values. When the self-part corrected scheme is applied 
to neutral hard spheres, using d^^ according to Eq. (32), 
also H^^{q) is predicted to very good accuracy even up 
to w 0.4, including the small-g region. 

We have shown that the (second-order) J7-scheme ex- 
pression for r]aa, given by Eqs. (25)-(27), is only applica- 
ble to hard-sphere like systems (i.e., at large screening). 
Its predictions for , however, agree well with the ASD 
data up to < 0.4. At larger 0, it underestimates the 
high-frequency viscosity. Based on arguments applying 
to low-salinity charge-stabilized systems only, we have in- 
troduced in Eq. (33) a simple correction to the (57-scheme 
result for rjoo ■ This corrected ^7 scheme predicts to good 
accuracy the high-frequency viscosity of low-salinity sys- 
tems, even up to the freezing concentration. Different 
from the self-part corrected S"f scheme for H(q) intro- 
duced in Eq. (31), which applies at any salt concentra- 
tion under the condition that (j) < 0.15, the corrected 
scheme for 7700 in Eq- (33) is valid at low salinities only. 
An appropriately corrected scheme for rjao , applicable 
for arbitrary salinities, is still missing. Its development 
could be the topic of a future study. 

An interesting application of the hybrid 5^ schemes 
for H{q) and rjao has been our validity tests of three 
approximate generalized Stokes-Einstein relations link- 
ing Tjao to ds, dcgei and dc, respectively. For an optimal 
test of these relations in the special case of neutral hard 
spheres, precise analytic expressions for d^^, H^^{qm), 
S^^{q„i), and S^^{q — >■ 0) have been used, valid in the 
whole fluid-state concentration range. The key finding 
from our validity tests of GSE relations is the strong de- 
pendence of their accuracies on the range and character 
of the particle interactions. The most striking example 
in case is the Kholodenko-Douglas GSE relation, which 
applies decently well to neutral spheres up to « 0.4. 
The very same relation, however, is strongly violated in 
low-salinity suspensions already at very low volume frac- 
tions. 
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